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ABSTRACT 

We investigate the contraction of accreting protoclusters using an extension of n-body 
techniques that incorporates the accretional growth of stars from the gaseous reservoir 
in which they are embedded. Following on from Monte Carlo studies by Davis et al., 
we target our experiments toward populous clusters likely to experience collisions as 
a result of accretion-driven contraction. We verify that in less extreme star forming 
environments, similar to Orion, the stellar density is low enough that collisions are 
unimportant, but that conditions suitable for stellar collisions are much more easily 
satisfied in large- n clusters, i.e. n ~ 30,000 (we argue, however, that the density of 
the Arches cluster is insufficient for us to expect stellar collisions to have occurred 
in the cluster's prior evolution). We find that the character of the collision process is 
not such that it is a route toward smoothly filling the top end of the mass spectrum. 
Instead, runaway growth of one or two extreme objects can occur within less than 1 
Myr after accretion is shut off, resulting in a few objects with masses several times 
the maximum reached by accretion. The rapid formation of these objects is due to 
not just the post-formation dynamical evolution of the clusters, but an interplay of 
dynamics and the accretional growth of the stars. We find that accretion-driven cluster 
shrinkage results in a distribution of gas and stars that offsets the disruptive effect of 
gas expulsion, and we propose that the process can lead to massive binaries and early 
mass segregation in star clusters. 

Key words: methods:Af-body simulations-starsdormation-stellar dynamics 



1 INTRODUCTION 

The formation of massive stars via collisions is an idea 
that has perhaps always been seen as exotic, compared 
to the less dramatic accretion model. In the context 
of clustered s t ar fo rmation, the idea was introduced by 
iBonnell et alj (|l998h . who considered the response of a 
star cluster as it accretes gas from its birth environ- 
ment. The idea was further explored in IBonnell fc Bate! 
(2002). However, the need for collisions as an alternative 
to accretion became less pressing as modelers of massive 
star formation pushed into two and three dimen sions (e.g. 
lYorke fc Sonnhaltedl2002l ; iKrumholz et~al]|2009t) , overcom- 
ing the radiation-pressure concern s of one-dimensional treat- 
ments (|Wolfire fc Cassinellilll987h . The success of accretion 
models, in addition to the extreme stellar densities required 
for collisions to become viable, pushed most consideration 
of collisions to the side. 

Early studies of collisional star formation tended to take 
the ONC as their touchstone to reality, a natural choice as 
the closest region of massive star formation. Orion by no 
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means represents an extreme of star formation; star clus- 
ters can be orders of magnitude more massive. More recent 
work has pointed toward these truly massive clusters as sites 
where accre tionaly induced col li sions may play a role in stel- 
lar growth. IClarke fc Bonne!! ()2008h developed analytical 
arguments, corrob orated by later Monte Carlo simulations 
l|Davis et al.ll20l5l ). suggesting that in populous clusters un- 
dergoing vigorous accretion, collisions in the core may be sig- 
nificant. Essentially, these arguments boil down to he fact 
that a larger-n cluster, with a longer two-body relaxation 
time, can be driven to higher densities by accretion before 
entering core collapse. 

Runaway collisions can occur in these clusters in the ab- 
sence of star formation considerations, driven by the gravita- 
tional dynamics of a system with a mass spectrum a ttempt- 
ing t o reach equilibrium (the Spitzer instability, ISpitzerl 
1969). Studies of this process i nclude both direct ra-body and 
Monte Carlo calculations (e.g. iPor tegics Zw art fc McMillan! 
120021 ; ICiirkan et all 12004 iFreitag et al.ll2006t ). A study tai- 
lored to the Arches cluster found that collisional runaway 
could take place on timescales of a few Myr, starting from 
a full mass functi on from 0.1-100 M(nin an initially mass- 
segregated state l|Chatteriee et all l2009h . It is very young 
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dense clusters like this that we are most interested in. In 
this paper, we work toward bridging the gap between hy- 
drodynamic studies of young gas-dominated clusters and 
purely dynamical studies of runaway collisions, using an re- 
body code extended to approximately treat the simultane- 
ous growth and dynamics of a massive embedded cluster. 
The space of cluster parameters and accretion scenarios is 
vast, and we do not attempt a comprehensive study, which 
will wait for follow-up papers. Instead we focus on idealiza- 
tions of two cluster types, one with modest numbers where 
collisions are not expected, and one more populous where 
collisions play a more prominent role. 



2 NUMERICAL METHOD AND INITIAL 
CONDITIONS 

The details of a star forming region can only be truly ex- 
plored with multi-physics, cluster-scale calculations. Cur- 
rently, the computational cost of these calculations precludes 
following the clusters out of the gas dynamical phase and 
into the purely stellar dynamical regime, and a simulation 
of an Arches-scale cluster seems well over the horizon. While 
n-body codes are much faster, cluster formation is an inter- 
play of gas- and stellar-dynamical processes, and the true 
applicability of purely gravitational studies in these matters 
is not clear. In this work we take a small step from the re- 
body side toward the multi-physics camp, by including a 
simple accretion effect in the cluster dynamics. 

The base cod e for our work is th e GPU-enabled ver- 
sion of NBODY6 l|Aarsethll2OO0l . I2003T I. We have modified 
the code so that stars grow in mass over some timescale; 
in this initial study at a constant rate proportional to their 
initial mass, so that the shape of the input mass function 
is preserved. In our current accretion prescription the gas 
that accretes onto a star is assumed to be at rest, so that to 
conserve momentum the stars slow down as they gain mass, 
and fall into the potential well of the cl uster. In contrast 
to the studies that mot ivated this work (|Clarke fe BonneTH 
120081 ; iDavis et aT]|2010t) the gas that accretes onto the sys- 
tem does not fall from outside the cluster. Instead, the initial 
cluster is dominated by a spatially-static potential identified 
with the natal gas. This potential is globally lowered at the 
same rate the stars grow, and is removed on an effectively 
instantaneous timescale at some point in the simulation, at 
which point accretion onto the stars halts and the experi- 
ment is a pure n-body one. 

At all times stars are assumed to have the main- 
sequence radius appropriate to their mass, which is likely a 
conservati vely small value (es pecially after a merger of mas- 
sive stars; ISuzuki et aill2007t ). When collisions are detected 
the stars merge with no mass loss. Winds and stellar evolu- 
tion are not included. The actual mass of the collision prod- 
ucts are thus a strict upper limit. The efficie ncy of a collision 
is realistically less than 10 per-cent (e.g. lLombardi et al.l 
ll996l ; lFreitag fc Bendl2005h . and the subsequent stellar evo- 
lution of a massive star can reduce the accumulated mass 
drasticall y, while increasing the radius (see the thorough dis- 
cussion in lGlebbeek et al"1l2009r ). A more realistic treatment 
of collisions would introduce further randomness to the sys- 
tem, which at this stage of these investigations we prefer to 
avoid. 
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Figure 1. Mass densities for run 32_2_a at four times, as labeled. 
Black histograms show the data from the experiment, with stars 
binned in equal- number bins. Gray curves show the value of the 
gas potential. 



We consider two types of clusters, focusing mainly on a 
populous cluster with n = 32768 (32k), as well as a lower- 
mass cluster with n = 2048 (2k). The details of the initial 
conditions are in table [T] The stellar initial mass function 
is a single power law with a Salpeter slope, £( m ) oc m -2,35 , 
over the mass range 0.03 - 3.0 Mq. Over 1 Myr, the stars 
gain mass so that they span the range 0.3 - 30.0 Mq, at 
which point the gas is expelled by fiat. All of the simula- 
tions presented here are performed with a global star for- 
mation efficiency of 30 percent. The stars and gas are both 
initially Plummer potentials with matching scale radii. Ini- 
tial experiments with alternative density structures do not 
qualitatively alter the results. These choices are arbitrary, 
although plausible; further papers will explore alternatives 
and enhancements to this basic setup. 



3 32K RESULTS 

Collisions are expected to occur in the approach to and af- 
termath of core collapse. In a single mass system with no 
gas potential, core collapse occurs via gravothermal con- 
traction. As energy is transferred from the interior of the 
cluster outwards via two-body relaxation, the stars in the 
center slow down, drop further into the potential, and gain 
energy. This negative heat capacity leads to a continuous 
rise in the central density as the c luster's outer radius ex- 
pands (|Lvnden-Bell fc Woodl ll968). As this is a relaxation- 
driven process, the ha lf-mass relaxa tion time is the natural 
timescale to consider (|Spitzerll 19871 ): 



t rh = 0.138 



1/2 3/2 

n ' r h < 



(Gm)V2in( 7 n) ' 



(1) 



where is the half-mass radius, n the number of stars, and 
m the mean stellar mass. The argument of the Coulomb 
logarithm depends on the mass spectrum of the stars, and is 
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Figure 2. Top row: Lagrangian radii for clusters with initial virial radii 0.75, 1.5, 2.5, and 3.25 pc, increasing from the left. Two stellar 
populations arc shown: in gray, the lightest third of the total mass; in black, the most massive third. The more compact clusters begin 
to mass segregate during the accretion phase, which ends at 1 Myr. Bottom row: the collisional history of the same clusters. The point 
at the apex of three points connected by lines is the merger product of of the other two points. Larger black points are those that exist 
at the end of the simulation, i.e. they are the end product of a cascade of mergers. 



usually taken in the range 7 =0.01-0.1, with the lower values 
for a ma ss-spectrum and the hi gher value for an equal-mass 
system (|Giersz fc Hegrid ll99^. Core collapse of Plummer 
sphere with equal masses takes place after about 15 t r h- 

When a cluster consists of stars with different masses, 
mass segregation drives the massive stars to the center. This 
massive core attempts to achie ve energy equ ipartition with 
the halo of lower-mass stars. ISpitzerl (| 1969T ) showed that 
for a two-component mass spectrum, the ability of the core 
of massive stars to achieve energy equipartition with the 
lower-mass halo depends on the relative masses and total 
mass of the two components. For realistic mass functions, 
the transfer of energy from the massive core to the halo ac- 
tually drives the system further out of equipartition, and 
core collapse ensues. This process, the Spitzer instability, 
generalizes to a mass spectrum instead of a two-component 
system (e.g. iGiirkan et al]|2004h . Core collapse with a broad 
mass spectrum generally occurs at t ^ Q.lt r h, much shorter 
than for the single-mass system. 

When the gas potential dominates the stellar potential, 
the radial scaling of the stellar cluster scales with the ac- 
cumulated mass as r oc M _1 , and the usual timescales for 
relaxation effects are complicated by the availability of a 
potential besides that of the collisional bodies. If the stellar 
potential begins to dominate a region, the accretion scenario 
is more akin to mass flowing onto a core from the outside, 
and the radius-mass scaling tu rns over to scale as r oc M~' i 
(these scalings are discussed in lBonnell et al.ll 19981 ). At this 
point familiar relaxation effects and the Spitzer instability 
are able to act, and the massive stars may decouple from the 
low mass stars and begin to go into core collapse. In figure 
[T] we show an example of the density evolution of the stars 
relative to the gas for one of our runs. Initially the gas dom- 
inates the stellar mass at all radii. As the stars gain mass 



and migrate inwards, the stellar density increases toward the 
center, eventually dominating the central potential. 

The behavior of our accreting clusters can be broadly 
divided into two regimes, depending on the initial condi- 
tions. In the first, the relaxation time of the core becomes 
short enough that the Spitzer instability can act prior to gas 
expulsion, so that both accretionaly driven contraction and 
dynamical effects drive the system rapidly to core collapse. 
In the second regime, the high and low mass stars remain 
coupled throughout accretion, which drives the contraction. 
When the gas is expelled the central regions are minimally 
affected, quickly revirialize, and the Spitzer instability alone 
completes the collapse. With our fixed mass accretion rate, 
the initial radius of the cluster determines which path the 
cluster evolution takes. 

In figure [5] we show the Lagrangian radii (Radii enclos- 
ing fixed mass fractions) for each of the four initial cluster 
radii that we examined. The two clusters on the left, 32_l_a 
and 32_2_a, are in the first regime where mass segregation 
has begun at the time of gas expulsion at 1 Myr. In the 
two larger-radii experiments, the mass segregation leading 
to core collapse commences only after the gas is expelled. 
We performed five realizations each at the cluster radii just 
on each side of the divide (series 32_2 and 32_3), to ensure 
that the qualitatively different collision behavior in these 
regimes is a general trend and not a quirk of a single run. 
We performed just a single run at each of the two extremes 
of our initial radius range, as we use them just to illustrate 
the trends in the bulk cluster properties, which vary little 
between realizations. The initial conditions and results are 
found in table [T] All the models were run for 3 Myr, ex- 
cept for run 32_l_a. This run became numerically problem- 
atic when the most massive star exceeded 600 Mq. By this 
point the simplifying assumptions we make regarding colli- 
sions have been stretched to their limits, and we stopped the 
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run at 1.6 Myr. The results quoted below are at 2 Myr after 
the simulation begins, to more clearly differentiate the im- 
portance of the accretion process from standard dynamical 
evolution. 



3.1 Collisional behavior 

The bottom row of figure [2] shows the collision trees for the 
clusters. In these plots, the product of a stellar merger is 
joined by lines to the two stars that took part in the col- 
lision. Collision products that are the end of an individual 
merger cascade, so are still around at the end of the simula- 
tion, are shown as larger black dots. In all but the largest- 
radius simulation, a clear runaway collision is seen as a a 
single star is involved in multiple, closely spaced collisions, 
rapidly gaining mass. In all cases, the seed for the collisional 
runaway is one of the most massive stars in the simulation, 
which is in a dynamically formed binary. The run 32_4_a 
the massive stars is formed from just 5 massive stars, which 
sequentially merge to form a 104 Mq object. However, this 
collisional growth only begins after ~ 2.25 Myr. 

While some collisions may occur prior to the minimum 
of the inner Lagrangian radii, in general the collisions begin 
in earnest once the core has collapsed and binaries have 
formed dynamically during three-body encounters. In all 
cases one or two stars rapidly accumulate mass to become 
the dominant body in the cluster, reaching masses in excess 
of 100 Mq. This process is similar to studies of runaway 
mergers leading to intermediate mass black holes mentioned 
in the introduction; the difference here is that this runaway 
may occur as a result of the cluster's initial growth rather 
than its later, purely dynamical evolution. The initial con- 
ditions of a standard n-body simulation would begin at our 
t = 1 Myr, when the more compact configurations have al- 
rea dy begun the collis i onal r unawa y process. 

IClarke fe Bonnelll <|200Sl ) and IPavis et al.l (|201(J ) pre- 
dicted that the number of collisions in an accreting cluster 
should be independent of the initial cluster radius, which is 
clearly not the case in this work. This difference is because 
those studies assumed that accretion onto the core contin- 
ues up through the point of core collapse. Because of our 
gas expulsion choice, this is only true for our run 32_l_a. 
An inspection of figure [5] shows that accretion results in ap- 
proximately an order of magnitude reduction in the cluster 
radius, and core collapse gives another order of magnitude 
for the inner Lagrangian radii, independent of the cluster 
radius at gas expulsion. The maximum density (and hence 
collision rate) in our models therefore depends on the initial 
radius. If we had allowed accretion to continue longer for 
the initially larger clusters, the relaxation time would have 
become short compared to the accretion time at approxi- 
mately the same density as the more compact clusters, and 
the peak density and collision rates would be similar. 

We re-emphasize here the idealized nature of the colli- 
sions in this work. There are a number of unmodeled pro- 
cesses that can alter the number of collisions and the final 
mass of the collision product. Our decision to take the radius 
of a merged star as its main sequence radius ignores the ini- 
tially out-of-equilibrium configuration, including a puffed- 
up envelope, that results from col lisions between main- 
seque nce stars(|Lombardi et alj|l996l ) or pre- main-sequence 
stars jLavcock fe Sills! 120051 ). An enlarged envelope can in- 



crease collision rates and the circularisation of hard binaries, 
and our collision rates are thus conservative. The most se- 
vere simplification we make, however, is n ot following the 
stellar evol ution of the colli s ion p roducts. Idebbeek et al.l 
(2009) and IChatteriee et alj |2009l ') showed, subject to un- 
certain modeling of stellar evolution at very high masses, 
that massive stellar winds can prevent the buildup of run- 
away products to more than about 100 Mq. The number 
of collisions per Myr in those authors' simulations were sev- 
eral times less than in these simulations, however, so those 
results can be taken only as a guide. 



3.2 Effect of binaries 

Dynamically hard binaries can serve to offset core collapse 
by inflating the cluster core in 3-body interactions. In that 
sense they work against a collisional scenario. However, a bi- 
nary can itself be perturbed sufficiently to lead to a merger 
of the two components. A full study of the effect of binaries 
requires a vastly expanded parameter space, but to illustrate 
the different behavior we might expect we repeat run 32_2_a 
with the 328 most massive stars (the most massive 1 per- 
cent) replaced by equal-mass binaries (run 32_2_a_b). The 
mass function of stellar systems is then the power-law IMF 
used in the other experiments, while the individual stellar 
mass function is different at the high end. This is not a re- 
alistic treatment, but serves to explore the general effects. 
The binaries are set up initially on circular, wide orbits; 
because the binary orbits shrink adiabatically just as the 
cluster does, we set them to be sub-AU scale at the time 
when core collapse begins, which places them on the hard 
side of the hard/soft binary divide by a factor of a few. 

The central density at 2 Myr in this run is virtually the 
same as in the equivalent run with no initial binaries. The 
'depth' of core collapse, measured by the minimum over all 
time of the inner Lagrangian radii, is deeper when binaries 
must be formed dynamically via three-body interactions. 
The cluster core inflates slightly in the latter case, while 
in the run with binaries the minimum is reached with no 
overshooting. The collisions in both cases are dominated by 
a single object, but in the run with primordial binaries more 
merger products are formed; 38 compared to 17 at 2 Myr. 
The binaries thus appear to act as the seeds of collision, but 
are no more effective at rapidly inflating the core than the 
dynamical binaries that form naturally. 

At the moment we do not include a treatment of tidal 
interactions, which can be an important source of binaries, 
although the efficiency of the tidal capture process is not 
fully understood. At densities of ~ 10 7 pc~ 3 , however, rate 
estimates for t hree-body capture (iGoodman et al.lll993T ) and 
tidal capture (|Press fe Teukolskvl 1977h indicate that the 
three-body channel will be more important in these cores. 
Another binary formation method is capture by a disc of 
material, which provides a dissipative environment around 
a star; the disc can result from tidal disruption of a low-mass 
star in a close encounter with a massive star (|Davies et al.l 
2006). This process is more important when a massive star- 
disc is involved, although the high velocity dispersion in th e 
core will decrease its effectiveness l|Moeckel fe Ballvl 120071 ). 
While these additional contributions to the binary formation 
rate could conceivably lead to an earlier onset of collisions, 
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the primordial population is likely to be the most important 
factor. 



3.3 Response to gas expulsion 

An interesting consequence of the shrinking scale radius of 
the stars relative to the gas is the robustness of the cluster 
to gas expulsion. By driving the stars to the center of the 
gas potential, the effective star formation efficiency is raised 
to levels such that the core of the cluster is virtually unaf- 
fected by the removal of the gas, even as the global efficiency 
remains l ow. A similar effect wa s seen (in a smaller-scale 
cluster) in lMoeckel fc Bate] l|2010h . where the long-term dy- 
nami c al evo lution based on the hydrodynamic simulation of 
iBatei (|2009h was less affected by gas removal than a simple 
estimate bas ed on the global star formati on efficiency would 
suggest fe.g. lBaumgardt fc Kroupall2007l ). This is an attrac- 
tive feature of this model, offering a way to begin the purely 
dynamical evolution in a very compact configuration despite 
the removal of the cluster's gas. 

The accretion scenario modeled here, with the gaseous 
background at rest relative to the cluster center of mass, 
yields the most extreme cluster shrinkage possible barring 
unrealistic possibilities such as a gas cloud with a counter- 
rotating stellar population. Examining figure [TJ we see that 
at just 0.3 Myr the inner regions of the cluster have reached 
an effective star formation efficiency of ~ 50 percent, and 
by 0.65 Myr the gas is negligible in the inner 0.25 pc. If the 
efficiency of the mass-loading of the stars were reduced, for 
example by locally-correlated gas and stellar motion, or the 
accretion phase represented a smaller fraction of the total 
mass growth, it still seems likely that the local star formation 
efficiency in the core will be quite high. 

The global star formation efficiency will change the de- 
tails of this effect, as well as the collisional behavior. How- 
ever, these changes should be small. When the gas potential 
dominates, the radial scaling of the stars is proportional to 
their mass as M~ , and thus the mass density of stars as 
M 4 . This is not changed by the amount of gas that makes 
up the potential. An inspection of figure Q] shows that even 
with 10 times more gas in the initial cluster, by 0.65 Myr 
the stars would begin to dominate the central potential. The 
time of core collapse would be pushed later, and since core 
collapse would begin from a more compact configuration, the 
peak density would be greater. As long as the stars reach the 
point of dominance, the resistance of the clusters to dissipa- 
tion from gas expulsion should be largely insensitive to the 
global efficiency. 



the density center, roughly the present-day tidal radius of 
the Arches. The peak densities in our models are of the 
order 10 6 -10 7 Mq pc -3 , higher than the estimated density 
of ^ 3 x 10 5 Mq pc -3 . This density estimate is based on 
the total mass of the cluster as extrapolated from counts 
of massive stars and a cluster radius given by the mean 
projected separation of each star from the centroid. We can 
form an analogous density estimate for our models, by taking 
the total mass and the average radius of each star. This 
value, as well as the observed estimate, are plotted as solid 
and dashed gray lines respectively. The agreement there is 
much better than comparing our peak density, although the 
uncertainty in the cluster mass relative to our model (Mtot ~ 
3 x 10 4 Mq) should be kept in mind. This mean density and 
radius are given as p m and r m in table [T] 

While the estimates of the Arches mass d ensity are 
limited to a mean value, lEspinoza et al. (2009) found an 
empirical surface density profile using stars in the range 
10 ^ A//Mq^C 120, using a King surface density profile, 

Their best fit for the Arches takes the values Eo = 2.2 x 
10 3 stars pc -2 , and r c = 0.14 pc. This fit contains ~ 535 
stars within 1 pc. While we cannot use the same mass range 
because of our arbitrary cutoff of 30 Mq , we can construct 
an analogous surface density profile with the most massive 
535 stars in our experiment. This is also shown in figure 
[3] The surface density of the massive stars in our model is 
clearly more centrally concentrated than the actual cluster. 

We are left with the following situation: in order to in- 
duce core collapse and collisional runaway before (or very 
shortly after) the accretion phase of the cluster, for reason- 
able accretion timescales the density of the stars at the age of 
the Arches is far too high to match observations. Cluster se- 
tups that more closely match the density profiles never begin 
mass segregation during the accretion phase; after the expul- 
sion of the gas, the cluster is similar to sta ndard initial con- 
ditions for a pure n-bod y simulation (e.g. IChatteriee et all 
120091 ; lHarfst et al] [2009). Accretion driven collapse as mod- 
eled here may have contributed to the compact observed 
conditions of the Arches, which could le ad to subsequent 
core c ollapse and collisions as argued by IChatteriee et all 
(2009). However, our simulations show that it is unlikely 
that the Arches was ever dense enough in the past to have 
already undergone a collisional runaway. The run that most 
closely matches the Arches, 32_4_a, still has a surface den- 
sity several times higher than observed, and it experienced 
only a single collision in the first 2 Myr. 



3.4 Density profiles: comparison to the Arches 

The present day Arches cluster is notable for its large mass 
and compact configuration. As the densest known cluster in 
the Galaxy, this is the natural place to look for the effects 
we are considering here. Mass estimates for the cluster range 
from ~ 1.2-2 x 10 4 M (|Figer et al.lll999l ; lEspinoza et~al] 
2009), and the average projected radius of the ob served stars 
is 0.19 pc (|Serabvn et al.lll998l ; linger et alj|l999h : this leads 
to a central density estimate of p c 3 x 10 s Mq pc -3 . 

In figure [3] we plot density profiles at 2 Myr for runs 
32f. For these plots we consider only stars within 1 pc of 



4 2K RESULTS 

Previous work suggests that collisions should be unim- 
portant at the number of stars and dens it y of Orion 
jBonnell et alj Il998l ; IClarke fc Bonnelll 120081 ; [Davis et all 
120101 ). As a check that our scenario does not result in sig- 
nificant collisions in such environments, and is able to yield 
a cluster with similar bulk properties to the ONC, we per- 
formed five realizations of a cluster with 2048 stars tailored 
to end up as an ONC-like cluster. The initial conditions and 
summary of the results for each of these runs are found in 
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Figure 3. Top row: volume density profiles at 2 Myr for the same runs shown in figure [2] The dashed horizontal line is the cluster 
density of the Arches found by iFiger et al.l l ll999h , and the solid horizontal line is the analogous density from our simu lations. Bottom 
row: s urface density of the 535 most massive stars in the simulation, for comparison to the empirical fit to the Arches from lEspinoza et al,l 
(2009), shown as a dashed curve. The gray histogram is at 1 Myr, the time of gas expulsion. The black histogram is at 2 Myr. 
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Figure 4. Lagrangian radii for the 2k cluster, shown as the me- 
dian value of the five runs. In gray, the lightest third of the total 
mass; in black, the most massive third 

table [T] The results quoted below and in the table are at 3 
Myr, roughly the age of the ONC. 

4.1 Collisional behavior 

In two of the five runs a single collision product is formed. In 
2b, a 14 Mq and a 6 Mq star merge, and this product later 
collides with the most massive star in the cluster. In 2e, the 
two most massive stars merge after forming a binary and the 
product later collides with a 0.3 Mq star. All of the collisions 
result from perturbed, dynamically-formed binaries in the 
very central regions of the cluster. 



4.2 Density evolution 

In figure 3] we show the Lagrangian radii of the 2k runs. 
Because of the noisiness of these data (due to the lower 
number of stars) we plot the median value of the five runs. At 
the time of gas expulsion the mass in the center of the cluster 
is dominated by the stars, just as for the 32k cases. The 
initial radius for these clusters is close to the 32_2 series, but 
because of the lower number of stars the relaxation time is a 
factor of ~ 2 smaller, and the runs are slightly further into 
Spitzer instability-driven core collapse before gas expulsion. 

The 2k clusters expands more after gas expulsion rela- 
tive to the 32k runs; nearly an order of magnitude expansion 
of the innermost Lagrangian radii, compared to a factor of 
~ 2 for the 32_2 series. This is not due to the removal of the 
gas, but rather the increased effect of a few hard binaries 
in the core when n is small. This r esult reinforces the ar- 
guments of lClarke fc Bonnelll ()2008l ) that stellar dynamical 
relaxation effects make the production and maintenance of 
high density cluster cores more difficult in low-rt systems. At 
3 Myr the central densities and half-mass radii of these mod- 
els compare w ell to the ONC's values of 2-3 x lO 4 Mq pc~ 3 
and ~0.8 pc l|Hillenbrand fc Hartmannl I1998T ). Just as for 
the 32k simulations the details of the final profile depend 
on the initial distribution, but this suffices to show that an 
ONC-like cluster can be created in this model without an 
alarming number of collisions. 

We also note that the result of the early core collapse 
in each case is a few binaries containing massive stars, and 
some degree of mass segregation (visible in the Lagrangian 
radii). It is possible that an accretion-driven core collapse 
process similar to the one modeled here could help resolve 
the mismatch between the apparent dynamical age of the 
ONC and its actual age. A n early dense phase in Or ion's 
history has been proposed bv lAllison et al. I (|2009| . |2010D as a 
way to form the Trapezium. Those authors invoke cold, frac- 
tal, gas-free initial conditions to drive the magnitude of the 
cluster's potential energy to high values, allowing collapse 
to a very compact configuration. The accretion scenario ex- 
plored here provides an alternative route to a similar result. 
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5 DISCUSSION 

The experiments presented here by no means constitute a 
complete examination of the importance of collisions dur- 
ing the accretion-driven contraction of forming clusters. The 
Plummer spheres used for the initial conditions and the 
highly idealized accretion model serve instead to illustrate 
the basic process. In future papers we will explore the more 
subtle aspects of this scenario, by adopting more realistic 
collision, accretion, and stellar evolution prescriptions, vary- 
ing the density profiles of the stars and the gas potential, and 
investigating the effect of binaries much more thoroughly. 

We find that the density profiles resulting from any clus- 
ter which hosts a significant number of collisions are too cen- 
trally concentrated to match known Galactic clusters. While 
the single- valued density estimate made in analogy to obser- 
vations is of the right order of magnitude compared to the 
Arches, the surface density is too high. The details of the 
density profile depend on the initial conditions, and a more 
sophisticated treatment of several processes including colli- 
sion behavior, stellar evolution, and the behavior of the gas 
in the simulations. Recall that the gas potential is fixed in 
radius and instantaneously removed. This choice was made 
for simplicity, and to avoid the temptation to massage the 
gas expulsion details to result in a particular density profile. 
Recently there has been some interest in fractal or otherwise 
clumpy initial con ditions, particularly in the context of early 
mass segregation (iMcMillan et al.ll2007i ; iMoeckel fc Bonnelll 
120091 ; I Allison et al.ll2009l , l2010l ). The effects seen in those pa- 
pers may accelerate the concentration of stars in the center 
of the cluster, potentially impacting the collision processes 
seen here. However since none of those works include any 
gas, which is almost totally dominant gravitationally dur- 
ing the formation stage we are exploring here, this is pure 
conjecture. In longer term ongoing simulations, we are ex- 
amining these factors in much greater detail. 

The basic results of these simulations seem to be ro- 
bust. In all large- n simulations that begin mass segregation 
prior to the end of accretion, many collisions occur, result- 
ing in a runaway object with a final mass well in excess of 
the maximum reached through accretion, 30 Mq in these 
models. In cases where mass segregation begins only after 
gas expulsion, the cluster evolution proceeds to core collapse 
as studied by previous authors. In these cases, their initial 
conditions correspond to the middle of our simulations. The 
tendency of one or two objects to dominate the collision 
process is important. This route is evidently not a means 
to fully populate the upper end of the mass function, but 
rather a means toward generating some of the most extreme 
objects. 

The similarity to previous work on the runaway growth 
of IMBH progenitors is clear. In this case, however, the core 
collapse that leads to runaway collisions is not necessarily a 
result of dynamics alone, starting from a cluster with a fully 
mature IMF. In our models, core collapse may be driven by 
the growth of the cluster. At the point when the natal gas is 
expelled the clusters are, depending on their initial radius, 
either already in core collapse or much closer to it than a 
typical equilibrium model of, for example, the Arches. As 
a result, collisions and runway commence on a very short 
timescale, usually within a few 10 yr. Importantly, the 
stars' characteristic radius shrinks relative to the gas po- 



tential, increasing the effective star formation efficiency and 
leaving the cluster resilient to disruption from gas disper- 
sal. Collisional matters aside, this presents a mechanism to 
create and maintain the dense initial conditions that a stan- 
dard n-body simulation of a dense young cluster might start 
from. 

Finally, we note that it is reassuring that this process, 
when tailored to reproduce a cluster with similar bulk prop- 
erties to the ONC, does not result in a significant number 
of collisions, despite going into core collapse prior to the ex- 
pulsion of gas. This result is due to the different behavior 
of low- and high-n systems when a few hard binaries form 
at core collapse. The fact that collisions are expected to be 
unimportant in Orion, the nearest and best studied site of 
massive star formation, is evidence against the need for colli- 
sions in order to create a typical O or B star. However Orion, 
much like current high-mass star formation simulations, has 
no single system in excess of ~ 50 Mq. It is the Universe's 
more extreme objects that may require more dramatic for- 
mation mechanisms. 
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Table 1. Parameters and results of the experiments. Results are at 2 Myr for the 32k runs, and 3 Myr for the 2k 
runs unless noted otherwise. 



Run name 


n 


ry.o/pc 


Pc /Mq pc~ 3 


rhm/pc 


Pm/M Q pC 3 


r m /pc 


Ncol 


M > 30 Mq 


oz_l_a (at l.D Myrj 


62 foo 


0.75 


5.1 x 10 


0.08 


2.0 X 1U 


0. 14 


57 


COO Ad on 

602, 46, 39 


32_2_a 


32768 


1.5 


1.2 x 10 7 


0.14 


8.0 x 10 5 


0.19 


33 


276, 37 


32_2_b 


32768 


1.5 


1.6 x 10 7 


0.13 


8.1 x 10 5 


0.20 


29 


344 


32_2_c 


32768 


1.5 


1.3 x 10 7 


0.14 


8.1 x 10 s 


0.20 


41 


356, 31 


32_2_d 


32768 


1.5 


1.4 x 10 7 


0.13 


8.0 x 10 s 


0.20 


26 


368 


32_2_e 


32768 


1.5 


1.3 x 10 7 


0.14 


7.9 x 10 s 


0.20 


26 


321 


32_2_a_b 


33096 


1.5 


1.3 x 10 7 


0.14 


8.3 x 10 s 


0.20 


80 


283, 46, 35 


32_3_a 


62 (68 


2.25 


O.O X 1U 


0.19 


4.2 X 1U 


0.24 


7 


100 


32_3_b 


62 (68 


2.25 


0.9 X 1U 


0.19 


A A w 1 a5 

4.U X 1U 


0.25 


10 


143 


32_3_c 


32768 


2.25 


6.0 x 10 6 


0.19 


4.2 x 10 5 


0.24 


11 


145 


32_3.d 


32768 


2.25 


6.6 x 10 6 


0.19 


4.3 x 10 5 


0.24 


6 


65, 45 


32.3.C 


32768 


2.25 


7.1 x 10 s 


0.19 


4.1 x 10 s 


0.24 


1 


57, 39 


32_4_a 


32768 


3.0 


2.5 x 10 6 


0.29 


2.1 x 10 5 


0.29 


1 




2a 


2048 


1.46 


4.5 x 10 3 


0.68 











2b 


2048 


1.46 


2.5 x 10 4 


0.65 






2 


50 


2c 


2048 


1.46 


1.7 x 10 4 


0.83 











2d 


2048 


1.46 


2.9 x 10 4 


0.58 











2e 


2048 


1.46 


2.4 x 10 4 


0.58 






2 


55 



N: number of stars in the simulation. 
rv,o'- initial virial radius. 
pc'. central density. 
rfmi- half-mass radius. 

Pm- for the 32k clusters, average density determined by the total mass and mean radius of stars within 1 pc. 
r m : for the 32k clusters, mean radius of stars within 1 pc. 
N co i: number of stellar collisions. 
M > 30 Mq: the masses of all stars above 30 Mq. 
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